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Abstract 

A brief review of the Fermi-Pasta-Ulam (FPU) paradox is given, together with its suggested 
resolutions and its relation to other physical problems. We focus on the ideas and concepts that 
have become the core of modern nonlinear mechanics, in their historical perspective. Starting 
from the first numerical results of FPU, both theoretical and numerical findings are discussed in 
close connection with the problems of ergodicity, integrability, chaos and stability of motion. New 
directions related to the Bose-Einstein condensation and quantum systems of interacting Bose- 
particles are also considered. 
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I. INTRODUCTION 



The goal of this paper is two- fold. First, we evaluate and summarize the most interesting 
results related to the FPU model, after the seminal paper |l| appeared in 1955. Second, 
we discuss new directions in the study of many-body chaos, that are related, directly or 
indirectly, to the FPU problem. We hope that our analysis will help future investigations of 
nonlinear classical and quantum systems of interacting particles. 

Numerous attempts to resolve the FPU paradox have resulted in a burst of analytical 
and numerical studies of nonlinear effects in physical systems. The primary interest of FPU 
was the observation of energy sharing in one-dimensional lattices with nonlinear coupling 
among rigid masses. For Fermi this study was directly related to one his first papers of 
1923 |2j in which he tried to rigorously prove the ergodicity hypothesis which lies at the core 
of traditional statistical mechanics. For a long time, the ergodicity was assumed to serve 
as the only mechanism needed for the foundation of statistical mechanics. Specifically, by 
assuming the ergodic motion of classical trajectories on the surface of constant energy, one 
can expect statistical behavior of a system and apply well developed statistical methods. 

In view of the result of Fermi [2j, and in accordance with wide-spread expectations, 
any weak nonlinear interaction between particles in a system with very many degrees of 
freedom causes ergodic behavior of a system. In fact, this point is used in the derivation of 
statistical distributions in the thermodynamic N —> oc limit for systems of noninteracting 
particles. Therefore, it is natural to expect that a system of 32 or 64 particles, as in 
the FPU numerical study, would reveal ergodic behavior, provided the nonlinearity is not 
extremely small. To the great surprise of FPU, the result was opposite: the long-time 
dynamics of the studied model appeared to be periodic, with almost perfect returns to the 
initial conditions. Having extraordinary intuition, Fermi noted that this effect may have 
very important consequences jj. About ten years later, two alternative explanations of the 
FPU paradox were suggested, giving rise to new phenomena, the integrability of nonlinear 
equations and dynamical (deterministic) chaos. 

One of the discoveries triggered by the FPU results was the complete integrability of a 
class of nonlinear differential equations. The first indication of this unexpected fact was due 
to numerical integration of the Korteweg-de-Vries equation in 1965 by Zabusky and Kruskal 
[J]. As was demonstrated, stable solitary waves (solitons) emerged from generic initial 
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conditions and traveled through the media, interacting with each other without losing their 
identities. This effect was suggested as an explanation of the remarkable recurrence in the 
FPU model, with the claim of its closeness to a completely integrable model. 

Another approach in resolving the FPU paradox was developed by Chirikov on the basis 
of his criterion of stochasticity (or dynamical chaos) p. As was already known from the 
study of nonlinear systems in applications to accelerator and plasma physics, the motion of 
a dynamical nonlinear system with few degrees of freedom can exhibit strong chaotic behav- 
ior. Thus, the description of these systems can be given in terms of conventional statistical 
mechanics. The mechanism of chaotic behavior of dynamical systems was found to be an 
exponential instability of motion for a wide range of initial conditions. The essential role 
in the emergence of this kind of instability is played by interacting nonlinear resonances. 
By 1965, the relatively simple criterion of resonance overlap enabled one to determine the 
conditions for the onset of stochasticity for various low-dimensional systems. Therefore, it 
was natural to apply the same approach to nonlinear systems of the FPU type. As a result, 
the threshold of stochasticity was found analytically in Ref.j^] and later confirmed numeri- 
cally 7}. According to these studies, the initial conditions used by FPU in their numerical 
simulations were chosen below the stochasticity threshold, just in the region corresponding 
to stable quasi-periodic motion. Above this threshold, the FPU model was shown to behave 
in accordance to the original expectations of FPU, revealing strong statistical properties, 
such as energy equipartition among the linear modes. 

In fact, the FPU study gave birth to a new method to study the physical laws of nature. 
The two first methods are well-known: theoretical and experimental physics. The new 
approach was predicted by Ulam and discussed in his mathematical book [8|. Expecting 
a future burst of computer technology, he proposed a new kind of synergetic cooperation 
between a physicist and a computer (see, also, Ref.j^]). Apart from the normal use of 
computers as a tool for the calculations of integrals, functions, differential equations, etc., the 
new role for computers consists of the study of physical systems ah initio, starting from given 
models and investigating their properties by varying parameters, forces, initial conditions, 
etc. This kind of activity was marked in the FPU paper as "numerical experiments." Later, 
this term has been widely used by Chirikov to stress the difference of the new approach from 
both theoretical and experimental studies. 

The FPU problem may be treated as a perfect example of this approach. Specifically, 
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first, the model was set up in the form of equations of motion. Then, "numerical experi- 
ments" were performed to see how rapidly the thermalization occurs due to the nonlinear 
interactions. And, unexpectedly, a new phenomenon was discovered, initiating further the- 
oretical studies. Afterwards, theoretical predictions were checked and further numerical 
studies gave new insight in the problem. Thus, the "synergetic" approach progresses. Since 
the time of the first "numerical experiments" of FPU many physical discoveries have been 
found first numerically, then explained theoretically and confirmed by real experiments. An 
exciting story of first twenty years of studies of the FPU paradox can be found in the book 
of Weissert 1101. 



II. THE FPU MODEL 

The primary aim of the authors of Ref.|j| was to observe thermal equilibrium in a non- 
linear string, and to establish the rate of approach to the equipartition of energy among 
different degrees of freedom. In order to treat this problem numerically, the continuum was 
represented by a large number of masses interacting with each other via nonlinear forces. 

The corresponding partial differential equations were approximated by a linear chain of 
particles of equal masses M connected by elastic springs. The linear part of the forces is 
determined by the constant K, resulting in harmonic frequencies loq = ^2K/M for all par- 
ticles. (Following many papers we assume, for simplicity, M = K = 1). For the nonlinear 
part, in Ref. |l| main attention was paid to the simplest cases of quadratic and cubic addi- 
tional terms, although some alternative forms of the interaction have also been discussed. 
For the quadratic force (called the a— model) the corresponding equations of motion are 

X n (^Cn+l 2x n -|- X n __x) -|- tt[(x n _j_i %n) (j^n x n—l) ]■ (-0 

Correspondingly, the chain of particles with additional cubic forces (called the f3— model) is 
governed by the equations 

x n = ( x n+l ~ 2a; n + x n—l) + P[{ x n+1 ~ x nf' ~ ( x n ~ x n-l) 3 ] (2) 

Here x n denotes the displacement of the n— th particle from its original position, and the 
parameters a and (3 are the strengths of nonlinear interactions between particles. 

In the absence of nonlinear forces the exact solution can be written in the form of normal 
modes Qkif) that are essentially the Fourier representation of the displacements x n (t) (for 
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fixed ends of the chain, xq = xn = 0), 

Qk{t) = ^ x n{t) sin ( 3 ) 

With this representation one can see that for any initial conditions x n (0) and x n (0) the 
energy 

E k = l -{Ql + ulQl) (4) 

of every k— th mode is constant. Therefore, the model is trivially integrable and energy 
equipartition among normal modes is not possible. As a result, the motion of this model is 
quasi-periodic in time, with the discrete spectrum determined by the normal frequencies u>k, 

UJk = 2sin {0)' ^ 

In the normal mode representation, the equations of motion take the form 

JV 

Q k + u 2 k Q k = a CijQiQj (6) 



for the a— model (JTJ) and 



Qk + <4Qk = P E DijiQiQjQi (7) 

i,3,l=l 



for the f3— model (J2J. Here Cij and -Dij,; are coefficients of the complicated dependence on 
the indexes i, j and /, defined by the nonlinear forces. 

Having in mind the predictions of conventional statistical mechanics, the authors of Ref. 
expected that by switching on the nonlinear terms in Eqs.Q-© [or equivalently Eqs.fJBJ- 
(0)], energy initially concentrated in a particular mode, will flow into all other modes, 
thus demonstrating the transition to equilibrium. In particular, analytical arguments were 
given in Ref. [ill] . according to which after a long time the systems of coupled anharmonic 
oscillators have to approach thermal equilibrium. Thus, it was a general belief that any 
kind of nonlinearity in a system with large number of degrees of freedom would give rise to 
ergodicity (see, e.g., Ref. [2]). And the latter was assumed to serve as the mechanism for the 
onset of statistical behavior in dynamical systems. 

The numerical studies in Ref. [lj were performed on Metropolis' new MANIAC computer 
with iV = 32 or N = 64 and with sufficiently small values of the nonlinear parameters a and 
j3, for zero initial conditions, Xq = x^ = 0. The first results refer to 1953-1954, with some 



additional runs that have been done after the death of Fermi in December 1954. Typically, 
the first mode with k — 1 was initially excited, for which most details are given. The time 
dependence of the energies E k (t) of all modes was studied for many fundamental periods 



T\ = 2-k/lui (for more details, see, e.g. Ref.|l2j). 

To the great surprise of the authors of Ref. Q], the results of a numerical simulation 
were quite astonishing. The behavior of both models was at first as expected: the energy 
spread to higher harmonics but after about 1000 oscillation periods T\, the flow of the 
energy into other modes stops, and the dynamics reversed, with the energy flowing back 
into the first mode. This recurrence of energy was found to be almost complete, with a 
decrease in energy of only about 2% of the total energy. In time, this periodic behavior 
persisted, thus demonstrating the absence of the expected statistical thermalization. The 
surprise was enhanced by the fact that the period of a recurrence was found to decrease 
with increasing coefficients of nonlinearity. Therefore, the nonlinear effects are significant 
and cannot be neglected. The time evolution did not lead to the equipartition of energy, 
rather, it demonstrated the existence of "quasi-modes" consisting of a number of linear 
modes. According to Tuck (j^j]), Fermi became really excited about this phenomena and 
thought that "something new and important might be at hand." 

A few possible reasons for this observed effect were discussed during the first stage of 
the story. Initially, the accuracy of the numerics was questioned, with a hint that more 
accurate calculations would show thermalization, although a very weak one. This point 
was somehow supported by the observation of a non-complete return of the energy in the 
originally excited mode. However, further more accurate computations of Tuck in 1961 (see 
Refs.0, revealed an even more exciting effect. It was found that at later times, the 
recurrence of the energy becomes more nearly complete. Specifically, a "super period" was 
found that is about 80 000 linear cycles, T\. The energy recurrence after this super period 
was more than 99% of the total energy. In g eneral, these results of Tuck, although not 
discussed in the literature until much later Jjj , confirmed the phenomena of the recurrence 
in the FPU-model. 

Of special interest was whether or not the energy recurrence can be associated with 
Poincare cycles that occur for ergodic systems. The estimate of the Poincare cycle for a 
chain of linear oscillators was derived in Ref. Q|. This estimate shows that the recurrence 
time in a chain of linear oscillators increases in an approximately exponential way with 
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the number of degrees of freedom. Therefore, it is clear that the Poincare cycles have no 
relation to the observed recurrence in the FPU-model. A careful analysis jl5[ of this estimate 
in application to the FPU-problem shows, however, that one should distinguish the FPU- 
recurrence from the Poincare cycles. The point is that the latter are defined for trajectories 
in phase space, rather than for the energy of a system. Obviously, the energy recurrence 
time will usually be much less than the Poincare time. Moreover, the estimate was given 
for a harmonic lattice, and there is no way, apart from direct computation, to determine 
this recurrence time in the presence of nonlinear coupling. This remark, however, does not 
change the conclusion that the FPU-recurrence has a different nature than the Poincare 
cycles. 

In view of the many discussions about the mechanism of irreversibility in the systems of 
interacting particles, it is instructive to mention some of the computations of Tuck. To see 
the influence of numerical errors, he performed the following check. After a few thousand 
cycles, the dynamics of the model was numerically reversed by the change of time and 
velocities of all particles. It was then found that 100% of the energy returns to the first 
mode. This fact was underestimated in the early 1960s. Now it can be treated as a direct 
(numerical) proof of the regular dynamics in the above models. 

As is now well known, dynamical chaos is characterized by an exponential sensitivity 
of the dynamics to the initial conditions. As a result, the unavoidable round-off errors 
in numerical simulations give rise to a drastic change for individual chaotic trajectories. 
Because of this exponential sensitivity and the round-off errors it is not possible to reach 
numerically the initial state, unless the reversal time is very small. This fact leads to the 
very important conclusion that chaotic systems cannot be treated as isolated ones since 
any weak external perturbation is essentially strong (see the discussion in Refs. [tgI Il7|). 
Therefore, this local instability serves as a mechanism for the apparent irreversibility in 
dynamical systems, although any dynamical system is reversible in principal (here, we do 
not discuss dissipative or noisy systems). 

For about a decade after the publication of the FPU preprint, discussions of the FPU 
paradox were restricted to a trivialization of the results, attempting to explain the recur- 
rence effect as simply as due to numerical errors, insufficient computation time, Poincare 
recurrence or the specific choice of nonlinear forces which prevents the ergodicity. It was 
still not well recognized that the FPU results initiated a new era in physics, associated with 
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nonlinear phenomena. 



III. PERTURB ATIVE APPROACHES 

The first analytical studies of the FPU paradox were described in the paper of Ford ^| . 
It was argued there that the absence of ergodicity in the FPU calculations may be due to the 
arithmetic properties of the unperturbed spectrum determined by Eq.(j5|). By making use 
of the perturbation theory of Kryloff and Bogoliuboff it was claimed that appreciable 
energy sharing among normal modes for a very weak coupling nonlinear interaction occurs 
only if the frequencies Uk of the unperturbed motion are linearly dependent (or, only if 
J2k m kUk = for some nonzero collection of integers {m^ = 0}). As for the FPU numerical 
data, they refer to the value of N as a power 2, therefore, to linearly independent frequencies 
(see details in Ref. 3)- For this reason, only few (low) modes in the FPU simulation could 
share the energy. Therefore, one should have multiple resonance conditions, in order to 
expect widespread energy sharing. However, as was shown in Ref.j^J, this idea, although 
quite useful in the description of weakly nonlinear oscillations, turned out to fail to explain 
the FPU paradox. Numerical experiments with many other values of N confirmed irrelevance 
of linear resonance conditions to the FPU recurrence. 

The numerical data of Ref.JlJ describe a relatively weak interaction between particles. 
This fact has triggered analytical studies of the FPU dynamics utilizing perturbation theo- 
ries. In an attempt to explain the quasi-period for the normal modes, in Ref. 2l( standard 
perturbation methods were examined in light of the application to long-time dynamics of 
nonlinearly coupled oscillators. As is known, the main problem is the small divisors that 
arise due to resonances between unperturbed oscillators. In the large N— limit the frequen- 
cies become dense and the frequency differences approach zero, causing all terms containing 
the small divisors to become infinite. Another problem is related to the appearance of 
secular terms that are proportional to a power of time t, and, therefore, restrict the appli- 
cation of time-dependent expressions to finite times. In order to avoid these secular terms, 
the Kryloff-Bogoliubov method [l^ was modified j^lj and applied j^cj to the FPU model. 
Specifically, second-order perturbation theory was found to give an accurate estimate (within 
15%) of the recurrence time and amount of energy exchange in the FPU problem. On the 
other hand, it was revealed that for some cases numerically studied in Ref. |l| , a higher-order 
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analysis is required due to a relatively large nonlinearity (when the nonlinear term is of the 
order of one-tenth of the linear term, in energy units). The important conclusion of these 
studies is that, strictly speaking, the FPU model does not belong to the category of weak 
coupling. 

It was also indicated jscj that when discussing the limit N — > oo, one should distinguish 
two different limits. The first one considered in Ref.j^ (see the discussion following), 
assumes that the length L = Na of the chain remains constant due to the decreasing 
spacings, a, between the particles. Correspondingly, the effective coupling a decreases with 
iV as 2a /N (correspondingly, (3 = 3(3 /N 2 ). In this way, by normalizing time t — > tN and 
the spacial coordinate z — > zL" 1 , one can obtain the following partial differential equations: 



and 



d 2 x 
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d 2 x 



dz 2 



(9) 



The corresponding initial conditions are prescribed over the range < z < 1 as x(z, 0) = 
Xq{z) and dx/dt\ t= o = 0. 

The other possible limit assumes the parameters of the chain are constant. Therefore, 
as iV — > oo, the length L becomes infinite and the frequency spectrum becomes dense. 
This is the limit typically discussed in the literature, especially, for the study of irreversible 
processes. One of the main questions is how the statistical properties of this system depends 
on the number N of interacting particles. 

As is shown in Ref . (stj , the phenomenon of recurrence in nonlinearly coupled oscillators is 
quite robust. Namely, asssuming the coupling constant in Eq.(JTJ) to be different for different 
particles (a kind of imperfection), one can "kill" the recurrence only with a sufficiently strong 
imperfection. This important fact indicates that the FPU recurrence is not an artifact of 
the chosen forces between particles. 

The meaning of the FPU recurrence and its relevance to the problem of ergodicity was 
thoroughly discussed in Ref. 



23| . Taking the results of FPU to be fundamental, it was 
suggested that ergodicity may not be required for the onset of thermalization in dynamical 
systems. In fact, this was a new insight into the problem of the foundations of statistical 
mechanics. In support of this point, it was noted that even a completely integrable system 
of linearly coupled oscillators shows a kind of thermalization for generic initial conditions. 



Specifically, after initially exciting one particular mass in the linear chain, one can observe 
an effective energy sharing among all particles. Having this analogy in mind, the authors 
of Ref . [23I performed an analytical study of energy sharing between linear modes in the 
a— model © with 2, 3, 5 and 15 oscillators. One of the results of this analysis was a 

n 

modification of the resonance condition obtained previously in Ref. 11 8| required for strong 
energy sharing, ^k^k ~ where are nonzero integers which depend on the particular 
coupling used. It was also found that, apart from this condition, strong energy sharing 
occurs for only certain initial conditions, not for all conditions. As a result, it was concluded 
that the dynamics of the FPU model may be consistent with the existence of additional 
integrals of motion; however, one can still speak about thermalization. To check this, the 
distribution of linear mode energies = Q 2 + io\Q 2 was obtained numerically for n = 5 by 
examining the time dependence of Ek{t). A very good correspondence to the exponential 
dependence was found, in accordance with the predictions of statistical mechanics. 

Thus, a new approach to thermal equilibrium problem was suggested: instead of the 
search of the ergodicity, one should study the conditions for strong energy sharing between 
normal modes. Moreover, it was suggested that in addition to the total energy, other 
integrals of motion may exist. At least one integral of motion was indicated to exist 23^ . 
due to a peculiarity of the model (the unperturbed part is purely linear, unlike many other 
examples for which in the absence of perturbation the unperturbed motion is nonlinear, see 
the discussion following). 

A new viewpoint according to which typical nonlinear systems are non-ergodic, has found 
rigorous confirmation based on the extensive mathematical studies of Kolmogorov, Arnold 

26], see, also, in Ref.j^). In 1954 Kolmogorov formulated 



and Moser (KAM theory, 



24, 



25 



the theorem that states that a weak nonlinear perturbation of an integrable system destroys 
the constants of motion only locally in the regions of resonances. In other regions of phase 
space, a set of points of positive measure remains for which quasi-periodic motion persists. 
This effect occurs for quite generic conditions on both the unperturbed motion and the 
type of perturbation. Loosely speaking, these conditions are as follows: the unperturbed 
system has to be nonlinear and the perturbation has to be weak enough and with a sufficient 
number of continuous derivatives. In the first proof of Arnold of Kolmogorov's theorem j^, 
for technical reasons, the number M of derivatives was assumed to be vary large. However, 
in subsequent studies by Moser [26], the minimal value of M was significantly reduced 
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(see, also, Refs. 



16 



17j | where corresponding estimates of this number were discussed and 



improved). Although a direct application of KAM theory to the FPU model is questionable 
(see the discussion following), the main result according to which one should not expect 
the ergodicity for a weak nonlinear perturbation, was essential for the acceptance of a non- 
ergodic dynamics in nonlinear lattices. 

Another perturbative approach has been suggested in Ref . j3| . It is based on the concept 
of Birkhoff-Gustavson normal forms as approximations to the Hamiltonians of the FPU 
lattices. Using these forms, one can show that for weak nonlinearity the motion of the FPU 
model is near-integrable. Further developments of this approach are reported in Refs. 29, 30] 
where the role of discrete symmetries and resonances was examined in great detail. In was 
also claimed that with the use of normal forms, the KAM theorem can be verified. 

In view of the KAM theory, it is important to mention the paper in which energy 
sharing and equilibrium were numerically studied for a chain of particles with elastic col- 
lisions. Specifically, in addition to linear forces between particles, elastic collisions were 
assumed caused by the finite diameters of particles. In contrast to the recurrence dynamics 
in the FPU model, in Ref.j^j] behavior that can reasonably be described as ergodic was 
found for N = 3 up to N = 32 particles. This ergodicity was ovserved both in the equipar- 
tition of the energy among all linear modes, in the time average, and by a rapid relaxation 
to equilibrium with the predicted values of temperature and pressure. This result was ex- 
tremely important for the understanding that the type of interaction between particles plays 
an essential role. As was understood later, the numerical data of Ref.|3l| are in agreement 

of ergodicity for a system of hard spheres with elastic 
collisions contained in a box. Moreover, apart from ergodicity, mixing was proved in Ref. 
which currently is considered as the most important property of dynamical chaos. Note that 
mixing automatically implies ergodicity. 



IV. INTEGRABILITY 



Inspired by the FPU paradox, in 1962 Zabusky analytically studied the continuous limit 
(JEJ) of the a— model He found an exact analytical solution for fixed initial conditions 
which turns out to break down at time t c ~ 1/e. At this time, the first derivative x z develops 
a discontinuity, therefore, x zz becomes singular. This means that the long-time dynamics 
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of the FPU cannot be approximated by the differential equations (JHJ). On the other hand, 
for time scales less than the breakdown time t c the comparison of the analytical solution 
with numerical data of Ref.[j| was reasonable. As was found in [22, 33 1, this critical time t c 
corresponds approximately to the time at which the energy in the second mode reaches its 



first maximum (if on 
Further analysis 



the mode with k = 1 is initially excited). 
331 ] showed that to study analytically what happens in the related 
continuous model, one must include the higher spatial derivatives that were omitted in 
taking the lowest continuum limit (JBJ)-©- To do this, one should use the following Taylor 
expansion of spatial and temporal differences: 



where a = L/N and x(z) = 
beta— model takes the form 



a 2 a 3 a 4 

zizQ,X z -\- ~~X ZZ i ~~~X ZZZ -\- %zzz • 
2 6 24 



(10) 



na. Therefore, the corresponding equation for the 



d 2 x 



d 2 x a 2 d 4 x 
Ih 2 + Udx 1 



(11) 



which descibes shallow water waves in classical hydrodynamics (see the discussion in 



Ref. 



For traveling waves in one direction only (e.g., to the right), one can approximately derive 
the equation 

Tr +u al + g W = ° (12) 

which is known as the Korteweg-de Vries (KdV) equation |35| . Here u = dx/d^, £ = 
z — cot, r = e*t, e* = \eco, and the velocity c in our units is 1. The parameter 5 2 = a/24/3 
is the dispersion which plays an essential role in discrete lattices. 

Apart from shallow water waves, the KdV equation is used to describe hydromagnetic 
waves in cold plasmas, ion-acoustic waves and long waves in anharmonic crystals (for ref- 
erences see Ref-lsj). Numerical study J^J of this equation in 1965 (see, also, Ref. 36]) has 
led to the discovery of solitary waves (or solitons), nonlinear waves that propagate through 
the media without changing their form. Specifically, it was observed that starting with the 
simplest initial condition [x(z, 0) = Csin7r,2 with x(z,0) = 0], solitons appear and strongly 
interact with each other, however, after interaction they preserve their identities. This re- 
markable property of stability of solitons (see, also, Ref. [37]) was treated as an indication 
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of the existence of a large number of integrals of motion. Due to the direct relevance of 
the KdV to the FPU model, the authors of Ref.j^J proposed that the recurrent dynamics in 
FPU lattices may be explained in terms of solitons as well. 

The numerical discovery of solitons |4j has attracted much attention to the KdV equation 

38| it was rigorously shown 



and triggered extensive analytical studies. In particular, in Ref. 
that two solitons keep their shape after interacting, and specific methods were proposed 
to prove the stability for the general case of any large finite number of solitons. Further 
heuristic methods were developed in Ref. {3^ to predict the number and speed of solitons 
emerging from arbitrary initial conditions. Later, a nonlinear transformation between the 
KdV equation and another nonlinear equation [namely, by changing the term uu z by u 2 u z in 



Eq. (jl2j) ] was found in Ref. 
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was developed in Refs. 
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After this, a rigorous method for solving the KdV equation 
42], by reducing the original nonlinear problem to a linear one. 
As a result, complete integrability of the KdV equation was proved for fixed boundary 
conditions. Similar properties of the KdV equation withperiodic boundary conditions have 
been found numerically in direct numerical simulations 36]. 

This remarkable integrability of the KdV equation was used to propose that similar 
properties may occur in nonlinear lattices. Thus, in Ref.j3] a nonlinear model (Toda- 
lattice) was introduced (see, also, Ref. [44]) with nearest neighbors interacting through the 



following potential: 



U(z) = %e- hz + az, (13) 
o 



where a and b satisfy ab > 0. The corresponding equations of motion have the form, 

(14) 



g -()(i„-i„-i) _ e -b(x n+ i-x n ) 



Formally, this lattice reduces to a harmonic lattice with the spring constant k = ab in 
the limit 6^0, keeping ab = const. One can also see that this model corresponds to the 
a— model when a = —b/2. Moreover, in the limit as b — > oo, the Toda lattice reduces to a 
hard sphere systems. Therefore, the model ()14|) covers two extreme limits of the interaction, 
from harmonic to hard-sphere. 

The first indication of the integrability of the Toda lattice appeared in numerical data 
of Ref. j^. It was shown that for N = 3 the trajectories cross the Poincare section in a 
way that the corresponding points lie on smooth curves. No evidence was found of regions 
with scattered points that is typical of integrable systems. Moreover, when studying the 
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divergence of neighboring phase trajectories, linear separations of the trajectories were al- 
ways observed. This fact is indicative of the stability of motion, unlike the opposite case of 
unstable motion for which the separation of trajectories increases exponentially with time. 
Similar behavior was found for N = 6 particles. Later, it was rigorously proved |46( that this 
lattice with periodic boundary conditions has N integrals of motion. Even the initial value 
problem can be solved for an infinite Toda-lattice by using the inverse scattering method, 
if the motion is restricted to a finite region of phase space; see details in Ref. 34j. It is 
important to stress that the integrability of motion in the Toda lattice does not prevent 
energy sharing among the linear modes (defined in the absence of nonlinearity) . As was 



shown in Ref. 



45| , energy sharing increases with an increase of nonlinearity. Therefore, good 



statistical properties can appear in completely integrable systems, provided the number of 
degrees of freedom is large, but not for all quantities. 

As a result of the very impressive discovery of integrability of the KdV and Toda lattices, 
it was often assumed (and this opinion still persists in some publications) that the FPU 
paradox was fully resolved by the concepts of integrability and solitons. However, the reality 
turned out to be even more exciting because of the direct relevance of the FPU problem to 
the dynamical chaos. 



V. STRONG CHAOS 
Analytical treatment 

Another approach to the FPU problem is based on the concept of dynamical chaos. 
For about ten years after 1955, an understanding of the fact that dynamical systems with 
few degrees of freedom may manifest quite strong irregular motion, turned into systematic 
studies of "stochasticity." This term was used to relate the irregular motion of completely 
deterministic systems to that known for physical systems which are governed by stochastic 
forces. Currently, two other terms are widely accepted, dynamical chaos, and deterministic 
chaos. These terms more correctly emphasize the purely deterministic nature of chaos. This 
is in contrast to conventional statistical mechanics which assumes, ad hoc, a probabilistic 
description of systems due to their underlying "randomness" . For a long time, the problem of 
establishing the conditions under which statistical mechanics is valid, was one of the central 
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problems in theoretical physics. The concept of dynamical chaos solves this problem, a fact 
which is still not well accepted by physicists although, is quite familiar to mathematicians. 
According to the modern viewpoint, classical statistical mechanics can be considered as a 
particular case of classical mechanics which describes both regular and chaotic dynamics. 
Thus, a statistical description, being useful and important, is an approximate one, and can 
be deduced from dynamical equations of motion under the conditions of dynamical chaos. 

One of the important studies of the problem of foundation of classical statistical mechanics 
is the work by Krylov j^. In his book, he analyzed the mechanism responsible for statistical 
behavior of dynamical systems, which is the exponential instability. By this term one means 
that the separation A(t) between two neighboring trajectories (in phase space) for generic 
initial conditions increases in time exponentially, A(t) ~ A(0) exp(ht). Here, the rate of 
the instability, h, is called the "dynamical entropy". Later, this quantity was rigorously 
studied by Kolmogorov and Sinai and it then assumed the name "Krylov-Kolmogorov-Sinai 
entropy" (KS- entropy). The positiveness of this quantity, h > 0, currently is used as the 
definition of dynamical chaos. Due to this instability, the dynamics of a system is extremely 
sensitive to its initial conditions, thus leading to mixing and other statistical properties (for 



details see, for example, Refs. 
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1Q: 

Numerical and analytical studies of dynamical chaos were strongly influenced by acceler- 
ator physics. As is known, the motion of charged particles in circular accelerators is affected 
by forces due to external magnetic fields that are required for the focusing of particles in a 
stable orbit. On the other hand, since the particles perform many revolutions (10 7 — 10 12 ) 
around a ring, nonlinear forces, although weak, are important for a long-term stability of 
particle motion. Early studies of nonlinear resonances in accelerators in 1956-1959 have led 
to the understanding that they can result in a kind of irregularity of motion. To the best 
of our knowledge, the first observation of this effect refers to the report 4||, in which the 
authors numerically studied the motion of electrons in a periodic electromagnetic field. Sim- 
ilar problems of stability have emerged when studying the motion of electrons in magnetic 
traps (see the references and discussions in Refs.|l6l Il7j]). 

The analysis of the stability of motion of particles in the presence of nonlinear pertur- 
bations has led Chirikov in 1959 to the concept of the overlap of nonlinear resonances j^. 
This term refers to the situation when nonlinear resonances strongly interact with each 
other. Specifically, it was found that when the nonlinearity is weak, one can consider any 
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particular resonance separately, making use of perturbation theories. However for a strong 

nonlinearity, the resonances cannot be treated separately because in the frequency space (or 

correspondingly, in action space) they are very close to each other. Thus, the overlap of 

resonances gives rise to the onset of a specific instability which leads to irregular (chaotic) 

motion. The analytic estimate for this overlap, known as the Chirikov criterion, turned 

out to be very effective for determining the conditions under which the dynamical chaos 

occurs in nonlinear systems ^?|. Results of the first experimental studies of the nonlinear 

resonances, their interaction and the onset of stochastic motion in electron-positron storage 

rings were reported in Refs. 

The application of the overlap criterion to the FPU model (J2J) has been reported in Ref . (f| . 

According to Eq.(J7J), there are two kinds of nonlinear terms. The term with i = j = I on 

the right-hand side plays a specific role and can be written separately, 

3/3 
AN 1 



Qk + ^\Qh 
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Here Uj km are the exact frequencies (including perturbation terms) of oscillations of normal 
modes that slowly depend on time. One can see that the selected term determines the 
nonlinear correction 

(5")k = -^l(2-4)<Ql> (16) 
to the linear frequency Even when small, this correction cannot be neglected since it 
depends on the energy of the k— th normal model. Note, that (5u)k can be obtained in 
first order perturbation theory, by averaging < Q 2 (t) > over the period of the unperturbed 
motion. 

The equations (|15|) describe the motion of nonlinear oscillators under the influence of 
external forces with amplitudes j3Fk m /8N. The spectrum of the perturbation due to these 

brces is given by the resonance frequencies u km — 2 sin n ^ k ^ m ^ with integers ±m = 0, 1, 2, ... 

f|. For small values of k <C iV (low modes corresponding to acoustic waves) the separation 
between resonances in frequency space is uj^+i — ^k ~ 2tt/N . Therefore, if nonlinear 
oscillations of a particular normal mode such as the maximal shift of frequency is much less 
than A w , then one can neglect the influence of neighboring resonances. In this case one can 
obtain the width, AQ, of a nonlinear resonance by keeping one resonance term only. The 
resonance criterion states that if AQ is of the order or larger than the separation 

A^ = u k+1 - u; k (17) 
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between neighboring resonances, then trajectories can no longer be associated with a partic- 
ular k— resonance and can wander between these two resonances. This transition from one 
to another resonance occurs in a quite irregular way, thus resulting in a kind of diffusion 
in frequency (or in action) space. In this case strong energy sharing between resonances is 
expected. 

For acoustic waves with A w m 2tt/N the critical perturbation for the resonance overlap 
is defined as ly], 



where E is the total energy of the lattice, therefore, E/N is the energy per normal mode, 
and Ak is the number of initially excited modes around the central k— mode. Here, units in 
which the distance between particles is fixed, a = 1, are used, therefore, the length of the 
chain is L = N. 

For the case in which high (optical) modes are initially excited, N — k <C N, the critical 
perturbation is given by the estimate 



Note that in this case the mean frequency spacing between nearest resonances is much 
less than that for low normal modes, A w « tt 2 /2N 2 , see Eq. f)17j) . For both acoustic and 
optical cases, the quantity 3/3 cr -j| ~ 3f3 cr (||) is the nonlinear term in the corresponding 
continuous model, which serves as a control parameter of perturbation. 

The above estimates determine the "stochasticity border" for the j3— model. From the 
condition (1181) one can see that for the lowest value k = 1 (the most studied case in the FPU 

ri 

model [1|), one must have a very strong perturbation in order to observe the non-recurrent 
behavior. Indeed, the nonlinear term in energy units in numerical studies of FPU has never 
exceeded 10%, therefore, the system was well below the stochasticity border. This explains 
the FPU paradox. On the other hand, from the estimate ()19|) for high modes with k close 
to N, the critical value of the nonlinear parameter f3 is relatively small and one can easily 
observe irregular motion with strong energy sharing among a large number of high modes. 
In fact, in a few runs of FPU with higher modes, one can see a more complex dynamics, and 
according to the expression (fT§|). the parameters used in Ref. [l| approximately correspond 
to the border of stochasticity. It was stressed in Ref.Q] that the border of stochasticity 
between quasi-periodic and stochastic motion is not sharp. Rather, it is relatively wide and 




(18) 




(19) 
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has a very complicated structure. For this reason, a strong dependence on initial conditions 
is expected in the transition zone. 

The analytical estimates obtained above refer to the overlap of two nearest resonances 
with k and k + 1. Therefore, they give the conditions for the onset of local stochasticity, and 
should not be treated as the threshold of widespread sharing of energy between all linear 
modes. The analysis [6] shows that when only a particular A;— mode is excited, the stochastic 
exchange of energy occurs to higher modes as well, at least for some group of initial linear 
modes. Indeed, with the flow of excitation in A;— space from low to higher k, the border of 
stochasticity decreases with increasing k, see Eq. (llSj) . It is interesting to note that, unlike 
the case of acoustic waves with k <C N, if high optic modes are initially excited, strong 
sharing between acoustic and optical modes occurs provided that low modes initially have 
a small amount of energy. 

Based on the resonance overlap approach, a specific study has been performed in Ref. 



for the a— model The analytical treatment has shown that this model appears to be 
much more stable than the f3— model. This is due to the fact that the nonlinear correction 
for linear frequencies in the first order of perturbation theory vanishes for the a-model [see 
Eq.©] and one needs to consider the second order approximation. Also, it was found that 
the most favorable conditions for the onset of stochasticity in this model occur when initially 
the modes with k ~ N/2 are excited. To compare with, both the limits of k -C N and k ~ iV 

n 

turn out to be more stable. In general, the analysis of Ref. |52| indicated that the a— model 
may be quite close to being integrable and further numerical studies have confirmed this. 

Numerical data 

Extensive numerical studies of the a— model have been performed in Refs.j^, 0, H ] 
exploring the above analytical predictions. As was 

property of stochastic motion is the exponential instability of trajectories with respect to 
a small change of initial conditions. To measure this instability, it was proposed to use 
the existence of the additional integral of motion in the FPU-model, namely, parity. As is 
clear from equations of motion, for fixed boundary conditions, Xq = Xn+i = 0, there is no 
interaction between the even, k = 2,4,6, and the odd, k = 1,3,5, modes. Therefore, 
when only odd modes are initially excited, the energy of the even modes has to be zero. 
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However, in numerical experiments 7} it was unexpectedly observed that when exciting the 
first mode with k — 1, the energy of even modes is not exactly zero, and moreover, this 
energy increases with time. Above the border of stochasticity the rate of this increase was 
found to be exponential, until the energy of the even modes approaches the energy of the 
odd modes. The mechanism of this phenomena is due to round-off errors (of the order of 
10~ 19 ) which cannot be avoided in numerical studies. One should note that these errors 
are not random since they are determined by a particular fixed algorithm. Therefore, they 
should be treated as a kind of dynamical perturbation which is not included in the original 
equations of motion. As a result, the energy of the even modes can be considered as a 
distance (in energy space) between two close trajectories. This concept was found to be 
very useful for determining the degree of instability. 

Given, after some initial time, a small amount of energy in the even modes (~ 10~ 14 of the 
total energy), the rate of instability for the j3— model (j2J) has been numerically computed as 
a function of the model parameters. Three regions of initial conditions have been examined: 
small modes with k = 1 or k -C N, high modes with k ^ N, and the intermediate region 
with k ~ N/2. In all of these cases quite good correspondence with the analytical estimates 
of Ref . (f| has been found. Specifically, for perturbations below the critical value, the rate of 
instability was approximately zero. This was in significant contrast to perturbations above 
the border, for which strong exponential instability was easily observed. In order to be sure 
that the numerical method used in determining the border is not an artifact, an additional 
check was done for two trajectories belonging to the same parity. The results were found to 
be analogous to those obtained from the energy increase of the even parity modes. 

In order to study quasi-periodic oscillations for normal modes involved in the dynamics, 



in Ref. 



55j | the possibility of describing a recurrence using truncated equations of motion 



was analyzed. It was found that for typical FPU conditions the dynamics of the model can 
be essentially described by a few equations for the modes close to the initially excited one. 
Namely, when exciting the mode with k = 15 for N = 32, three (with k = 14 — 16) and five 
( with k = 13 — 17) coupled equations have been examined numerically. 

The important question studied numerically in Ref. Q is the dependence of the rate h of 
instability on the perturbation parameter (3. As expected from the predictions of Ref. 
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for large values of (3 the dependence 



/i^A^ln-^ (20) 

Per 



has been found to correspond to numerical data, with A w as the mean distance between 
the unperturbed frequencies, see Eq. ([17j) . On the other hand, when the perturbation was 
not very strong, the dependence turned out to be very different and it can be fitted as 
h ~ A w (/?//3 cr .) 4//3 . It was proposed that for a weak enough perturbation, the instability, 
being exponential, is due to high order resonances. Although these resonances are more 
dense, the diffusion among these resonance is much slower. Therefore, apart from the strong 
stochasticity (chaos) determined by the overlap of main resonances (due to the first order of 
perturbation theory), one can speak about weak chaos which can also lead (on much larger 
time scales) to strong energy sharing between normal modes (see Sect. VII). 

In addition to the instability of motion, in Ref.Q other statistical characteristics of the 
dynamics have been studied as well: energy sharing among modes, the time dependence of 
the energies of each mode, time correlations < x n (t)x n (t + r) > and < Ej,(t)E^(t + r) > for 
displacements and energies, as well as correlations between energies of different modes. The 
results strongly support the onset of strong chaos above the border as analytically predicted 
in Ref.fl. 

An additional numerical study has been reported in Ref . [5^ (see, also, Ref.jsjj]) with 
higher accuracy and a larger number of particles (up to N — 500). These new data confirmed 



the main findings of Ref. |7J concerning the border of stochasticity. Moreover, the exponential 
dependence (J2*U|) for the rate of exponential instability was also supported by these data. 
One of the new observations was the existence of an initial time scale on which a non-chaotic 
excitation of modes different from the initially excited one occurs. After this initial time, a 
stochastic exchange of energy begins. Therefore, it was proposed to modify expression (fTHj) 
for the stochasticity border by using larger values of k due to this effect. This effect seems 
to be relevant to the emergence of solitons that occurs in the Toda lattice. Therefore, the 
following picture seems to be more correct: first, an initial regular dynamics occurs in the 
model, with the excitation of higher modes. After some initial period of time, a stochastic 
exchange between the modes comes into play, with a practical irreversibility of motion and 
onset of thermalization. Note that this thermalization can be restricted to a finite number 
of modes, much less than the total number N of degrees of freedom. 
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To compare with the KdV equation (JT2J), one should recall that this equation is an 
approximation to the FPU model. The main difference lies in the assumption that the 
waves traveling in different directions can be considered independently. This fact may be 
crucial in the analysis of the application of the KdV to the FPU model. In order to check 
how important the above approximation is, in 

Refs.0 Q a specific study of the FPU 
model with periodic boundary conditions xq = xn was carried out. The quantity of interest 
was the emergence of waves traveling in a direction opposite to that of the initial wave. 
Specifically, the percentage of energy of standing waves in comparison with the total energy 
for the k— th mode was calculated as a function of the perturbation f3 at some (large) fixed 
time. It was found that for small perturbations, opposite moving waves practically do not 
appear. However, for (3 ~ f3 cr the amount of energy in the waves in the opposite direction 
is of the order of total energy. This fact demonstrates that for large perturbations the FPU 
model is very different from the KdV model. 



VI. FURTHER RESULTS 



Energy sharing and equipartition 

The existence of an initial period of time ( "induction period" ) for which the motion does 
not reveal strong energy sharing, was studied in detail in Ref. J3] for the (3— model with zero 
boundary conditions [see, also, the study [58] of the 2D model]. After this period, strong 
energy sharing between a large number of modes was clearly seen, thus corresponding to 
the predictions of Refs.jO] of the onset of strong stochasticity. It was also found that this 
period increases as the nonlinear coupling decreases. As the criterion for the establishment 
of thermal equilibrium, in Ref.j^Z] velocity- velocity correlations between close in the chain 
particles were used. It was shown that below a critical value of the nonlinear coupling these 
correlations are very small; this was used as an indication of thermal equilibrium. These 
results were compared with those obtained for the model with linear coupling only (k| for 
which the correlations were found to be stronger due to the absence of ergodicity. 

However, the approach of Ref. ^3 based on the examination of correlation functions has 
been criticized in Ref. 0] (see, also, Ref. (§]]). It was argued that this method does not give 
global information about the phase space of the system, and strongly depends on the choice 
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of correlation functions. Analyzing other methods for determining the chaotic transition, the 
authors of Ref . j^] proposed to study the distribution of energy modes after a relatively short 
period of time. Their analytical analysis showed that at short times this distribution has an 
exponential dependence, W(k,t) ~ exp[— B(t)k], on the function B(t) that depends on the 
model parameters. Numerical data have shown that with a high accuracy the distribution 
of energies corresponds to the analytical expression for B(t). At later times for large enough 
nonlinearity the distribution W(k,t) was found to be of the expected form W(k,t) ~ 
corresponding to the Boltzmann distribution of mode energies. These results also confirm 
the existence of the stochasticity threshold in its dependence on the nonlinearity parameter 

An interesting quantity to measure the energy sharing has been proposed in Ref. [62]. 
This quantity was found to be quite useful in the study of relaxation properties of nonlinear 
lattices. In order to characterize the energy spread in the mode representation, the spectral 
(Shannon) entropy S(t) is used, 

TV 

S(t) = -£>*(*) In «;*(*), (21) 
fc=i 

where Wk = Ek/J2iEi is the normalized energy of a particular mode. The spectral energy 
is zero when only one normal mode is excited, and reaches its maximal value S max = 
ln(iV/2) for complete equipartition of the total energy among all modes. To avoid the clear 
dependence on the number of oscillators, the normalized quantity 

^ Smax (22) 
Smax ~~ S(0) 

was introduced. Here Soo is the maximum value reached by S(t) in the time evolution, 
associated with the "asymptotic" value of S(t). This quantity r\ is bounded between zero 
(which corresponds to complete "localization" in one mode), and one (which corresponds to 
perfect equipartition). 

Computing the normalized spectral entropy rj, in Refs. (3, 1^ strong evidence in favor 
of the existence of an equipartition threshold was given. In numerical simulations, periodic 
conditions were used for the (5— model, with initial excitation of a group of modes k ± Ak/2 
with small values of k <C N. The integration period was chosen large enough to ensure a 
practical independence of r\ on time. These numerical data revealed the remarkable result 
of a universal form of dependence of r\ on the energy density e = E/N for many values of 
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N from 64 to 512. It is interesting to note that this effect is insensitive to randomization of 
the FPU model, for which linear forces were assumed to be different for different particles 
(see details in Ref.jo^). According to these results, the threshold of equipartition does not 
disappear in the large N— limit. 



Since in numerical computations 



62 



63| the mean value of k for initially excited modes 
was taken to be proportional to iV (as well as Ak ~ N ), it was claimed, that the ob- 
tained results are in formal contradiction to the condition (118)1 of a strong chaos, where the 
threshold disappears with Ak ~ k ~ iV — > oo. The explanation of this contradiction lies 
in understanding the meaning of the stochasticity threshold (fTHJ) . Indeed, according to its 
derivation, this condition refers to the overlap of nearest resonances only, and there is no 
direct relation to energy sharing among all modes. Although the numerical data show that 
a large number of modes appears to share their energy above this threshold, the question 
about complete equipartition remains open. 

An important comparison of the FPU model with the Toda lattice is given in Ref. 64]. As 
was already discussed (see, e.g., 

Refs.HQ), 

quite strong energy sharing may be observed 
in the FPU model in the regime of strong recurrence, below the border of strong stochasticity. 
For this reason, the dynamics of the Toda lattice and FPU models were analyzed from the 
viewpoint of equipartition. As was pointed out, one should distinguish between "energy 
sharing" and "equipartition". It was shown that strong energy sharing can be observed in 
both models. However, equipartition occurs in the FPU model only, which is understood to 
be non-integrable. The numerical data which demonstrate this difference are based on the 
form of the normalized spectral entropy r\ for large times. Namely, for the Toda lattice the 
spectral entropy Sit) never reaches its maximum value, in contrast to the FPU model for 
which it does reach the maximum. Therefore, for the Toda-lattice the dependence of r\ on 
the energy density e = E/N does not show a transition to zero. 



Stability conditions 

In order to characterize the difference between integrable and non-integrable lattices, in 
Ref. 64] it was proposed to study these models from the viewpoint of stability. Specifically, 
it was observed that the time dependence of the trajectory in the artificial phase space 
rj(t), r](t) is clearly different for these two cases. Were the FPU trajectories to appear 
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irregular and unstable to an external perturbation, for the Toda lattice the trajectories 
reveal clear quasi-periodicity and stability to this perturbation. Thus, in order to distinguish 
between integrable and non-integrable chains, a study of the stability of motion is needed. 
The first analytical study of the stability of motion in the FPU model is reported in 



Ref. 55]. It was found that the recurrence dynamics in the (3— model with fixed boundary 
conditions may be correctly described by keeping a small number of equations for those 
modes that are essentially involved in the dynamics. For these equations, one can write 
the condition of a linear stability which stems from the corresponding Mathieu equation. 
According to this condition, a critical value of perturbation exists above which the motion 
is unstable. However, the relevance of this instability to the overlap of nonlinear resonances 
remains unclear. 



Another approach to instability of motion has been developed in Ref. [66j in application 
to the KdV equation. It was analytically observed that certain KdV solutions are unstable. 
Using this fact, an attempt was made to relate this instability to that observed in the FPU 
lattice. The analytical predictions obtained for the KdV model have been claimed to explain 
the instability of motion in the FPU model. In this study specific initial conditions in the 
form of a cnoidal wave were used, for which numerical data manifested a good correspondence 
to analytical predictions. 

The above analytical studies have suffered from the absence of reliable expressions that 
would predict the dependence on the nonlinear parameter and the number of particles. The 
first attempt to shed light to this problem was made in Ref. where the stability condition 
was obtained for the specific case of the highest linear frequency, which is initially excited. 
Note that for periodic boundary conditions the linear spectrum is doubly degenerate. There- 
fore, the highest frequency corresponds to the middle of the spectrum, k = N/2. Using a 
variational equation for this mode, a simple approximate formula for the f3— model (in the 
large N— limit) was derived, 

3 A | - H. (23) 

Applying this estimate in Ref.|57( to numerical data obtained for iV = 15 and zero boundary 
conditions, some discrepancy was found. In this respect the authors claimed that their 
critical value is, in essence, an upper estimate for the threshold of the chaotic transition, and 
cannot be compared with the condition of widespread energy sharing. Another uncertainty 
is due to the different boundary conditions used in the analytical evaluation. Later, the 
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stability condition (|23|) was obtained [68| (with almost the same constant) in the more 
general context of bifurcations of periodic orbits in nonlinear Hamiltonian lattices, and with 
some relation to symmetry breaking and vibrational localization (breathers). 

A similar analysis has been done for the a— model. The corresponding stability 
condition turned out to have the same form, 

E °- 84 

a -N " -W- (24) 

This result is quite unexpected since the a— model is assumed to be closer to being integrable 
than the (3— model. An additional study of the stability of the a and (3— models with 



attractive potentials, a < and f3 < was performed in Ref. 67]. 



The important results are reported in Ref. 



691 ] . Using the so-called narrow packet approx- 



imation, the following stability condition was derived for the f3— model: 

%f « w (25) 

It was found that for f3 > (3 g a parametric instability of motion emerges for wave packets that 
populate a number of linear modes for k within the interval \k — k \ = Ak -C ko in a region 
of the optical phonon spectrum around k « N/2. Note that for the periodic boundary 
conditions used in Ref. [69] . the value k = N/2 corresponds to the mode with highest 
frequency, due to the double degeneracy of unperturbed frequencies u k . The point is that 
for such initial conditions the FPU model reduces to the nonlinear Schrodinger equation, 
which is known to be completely integrable (see details and discussion following). This fact 
establishes a link between nonlinear lattices of the FPU type and models which are now 
widely used in application to the Bose-Einstein condensation. 

According to the condition (J25J), if the amplitude of the initially excited modes exceeds 
some critical value, the parametric instability results in a rapid spread of the wave packet 
over many linear modes. Even though this packet spreads rapidly, the narrow packet approx- 
imation remains valid for some time. One can see that formally Eq. (j23j) coincides with the 
condition for the onset of stochasticity due to the overlap of two close nonlinear resonances, 
see Eq. (jl9j) . It should be stressed that, strictly speaking, both conditions correspond to the 
low boundary for the emergence of chaos and may be different from those for strong energy 
sharing among all modes. Note also that the result of Ref.|]§9f practically coincides with 
Eq.([23|h which also refers to the instability of the highest mode. 
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It is also very instructive that the condition of the validity of the narrow packet approx- 
imation obtained in Ref.j3| reads as 

E 7T 2 

*>*S " N- (26) 

This condition may be treated as the critical value above which strong equipartition among 
all modes arises in the lattice. The important point is that the additional factor iV stands in 
the denominator of Eq. 1)26)1 in comparison with the stability condition l)25j) . Therefore, one 
can propose that the difference between the critical value of perturbation for local (chaotic) 
energy exchange between nearest modes in the k— space, and global equipartition of energy 
in the lattice, is mainly due to this additional iV— dependence. Note, however, that what 
we discuss here refers to the initial excitation of the high frequencies only. The problem of 
stability for small values of k (acoustic waves) seems to be very different. 

The problem of stability of solutions in the /3— model with periodic boundary conditions 
has been studied in a generalized approach in Ref.jzi|. A complete rigorous analysis has 
been done for the 7r— mode considered in Ref.j^}, resulting in an additional correction term 
which depends on the number of particles. It was shown that other exact solutions exist 
below some critical value of nonlinear parameter (3. Moreover, the presence of multi-mode 
invariant manifolds was shown. It was also pointed out that the relevance of the instability 
of these solutions to widespread stochasticity is not clear, although numerical data generally 
manifest such a connection [compare Eq. (J23j) and Eq. (j25J) with Eq. ()19| where k « N). The 
general case of a nonlinear lattice with both a and (3 terms has been under careful study in 
Ref . [li | where, in particular, the relevance of the instability of the highest frequency mode 
to stable localized solutions ( breathers) has been discussed (see, also, Refs.|72l. l73j|). The 
role of periodicity of boundary conditions can be examined with the use of Birkhoff normal 
forms, see details in Ref.jz^j]. 



Lyapunov exponents 

As is discussed in Ref. 0, the important quantity that characterizes dynamical chaos, is 
the local instability for which two close trajectories in phase space diverge, in time, expo- 
nentially fast. For this reason, many modern numerical studies are based on the calculation 
of the rate of this instability (KS-entropy). A detailed analytical and numerical analysis 
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of the method of calculating the KS-entropy in many- dimensional dynamical systems has 
been performed in Ref. The approach developed by these authors was found to be 
extremely useful in the study of dynamical chaos in various physical systems. Extensive 
numerical analysis of KS-entropy in application to anharmonic chains with a Lennard- Jones 
interaction is reported in Ref. 76]. This lattice is known to exhibit the stochasticity tran- 
sition in a more clear way, compared with the FPU model (see, e.g., Ref. 77[). These data 
show that the standard procedure of computing the dynamical entropy due to the average 
h ~< In \d(t) / d(0)\ > along the trajectories is quite stable with respect to different kinds of 
computational errors, and gives reliable results. However, it should be noticed that in this 
method the quantity h is not exactly the Kolmogorov entropy, although it often is close to 
it. The exact expression for the KS-entropy is a sum of all positive Lyapunov exponents 
(LE) (for details and references, sec, e.g., Ref. 48]). Actually, the above method determines 
the largest Lyapunov exponent, and this is sufficient to distinguish between chaotic motion, 
h > 0, and (quasi)-periodic, h — 0, motion. On the other hand, it was argued in Ref. 0] 
that the largest Lyapunov exponent may not show a correct picture since it does not ex- 
perience different ergodic regions. However, in Ref. [3] an opposite conclusion was drawn 
from numerical data, according to which different non-connected regions of chaos can also 
be detected, by searching the fluctuations of the largest Lyapunov exponent. 

Since the largest Lyapunov exponent Ai ~ h can be used as a measure of chaoticity in 
a system, analytical estimates of Ai and its scaling properties are extremely important. In 
Ref . it was argued that the spectrum of Lyapunov exponents for long chains may be 
well approximated by the Lyapunov exponents of products of independent random matri- 
ces, provided the energy per mode, e, is sufficiently large. This point has been used in 
Ref. 80] where an analytical estimate for Ai was derived in the large A^— limit. Specifically, 
a Gaussian model with noise was used as an approximation of the dynamical FPU model, 
and the analytical results were compared with numerical calculations. An amazingly good 
correspondence was found between analytical predictions and numerical data, and two scal- 
ing dependencies were dicovered, Ai(e) ~ e 2 for e — > and Ai (e) ~ e 1 / 4 for e — > oo. In 
comparison with the previously obtained numerical data in Ref. [81 1, the first scaling was 
confirmed. As for the second one, for large e, their result Ai ~ e 2 ^ 3 of Ref.jsiJ was different. 
A later theoretical study in Ref. 82| confirmed the dependence At (e) ~ e 1 / 4 by making use 
of a different approach. An important point of the studies in Refs. [sol Isil . Is^ is that there 
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is a clear transition from one scaling to another. According to Ref. 821 ]. this transition oc- 
curs at E pd n 2 /3N/3 which corresponds to the onset of strong chaos, see (jl9j) . This very 
important fact confirms previous findings that below the border of strong chaos, found from 
the Chirikov criterion, weak chaos persists. Thus, one may expect energy equipartition for 
any weak nonlinearity, however, after much longer times (see discussion below). 

Much more information can be drawn from the knowledge of all Lyapunov exponents, not 



only from the largest one. A numerica 
systems has been developed in Ref.js" 



method of computing all LE in many- dimensional 



8J]. Currently, this method is widely used in many 



applications, both classical and quantum. Of particular interest is the distribution of LE 
as a function of the index j according to which the Lyapunov exponents are ordered in 
an increasing way. Spectrum of the LE has been numerically studied for the (3— model in 
Ref. 3] where it was found that already for 40 to 60 particles the limiting distribution 
emerges. Thus, these results demonstrate the existence of a thermodynamical limit for the 
spectrum of LE. The obtained distribution was discussed in Ref. 185] with a view towards its 
relevance to random matrix approaches. 

A very interesting study was reported in Ref. [sfj. As was already mentioned above, the 
a— model seems to be more stable than the f3— model. For the first time this point was 



noted in Ref. 



521 ] where the resonance overlap criterion was obtained for the a— model. The 



same conclusion can be drawn from the analysis of nonlinear differential equations of the 
KdV types, from the viewpoint of their integrability 0]. In order to quantify the difference 
between the integrable Toda lattice and non-integrable a— model, in Ref. 86] the largest 
LE, Ai, was computed for both models as a function of time. It was found that at large time 
scales which are inversely proportional to the energy density e = E/N, the time dependence 
of Ai(t) is practically indistinguishable for both models. However, starting from some critical 
time, a drastic difference is clearly seen: for the a— model Ai tends to a positive value, in 
contrast to the Toda-lattice where Ai continues to vanish with an increase of time. As a 
possible explanation of this remarkable effect, the authors conjectured the coexistence of 
tiny chaotic regions and relatively large regions of stable motion in the phase space of the 
(3— model. It was argued that the trajectory might be trapped for a long time in a relatively 
large stable regions. Then after some time, the trajectory "finds" a way to leave the stable 
region and enter a stochastic region. Apart from this suggestion, there is no satisfactory 
explanation of the observed effect. Note that the computation performed in Ref. js^ was in 
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exact correspondence with the old computations of FPU, however, with the highest accuracy 
and for the longest time ever achieved. Specifically, the lowest frequency mode with k = 1 
was initially excited. 

The observed effect is interesting from many viewpoints. First, it again confirms the 
point that the practical difference between integrable and non-integrable models may be 
very small, and this difference can be detected only on very large time scales. Second, it 
gives direct evidence of the existence of the threshold of weak chaos, although indirect 
indications for the (3— model have been reported before (see, e.g. Ref.|88(). Third, it shows 
that the largest Lyapunov exponent Ai seems to be the only quantity which can give reliable 
results concerning weak chaos. By searching other initial conditions, the authors of Ref. 8(| 
claim that the existence of the stochasticity threshold e c can be clearly seen by examining 

Another important result of Ref. [8o| is the N— dependence of the stochasticity threshold 
for the case when all modes are initially excited. With an increasing number of particles, 
it was found that the scaling dependence e c ~ 1/N 2 is well supported by the numerical 
data. Thus, in the limit N — > oo the threshold in the FPU model vanishes. So far, there 
is no satisfactory explanation of this result. Note that most studies have been done for the 
(3— model, and it is questionable whether there are quantitative properties shared by these 
two models. 



VII. WEAK CHAOS 

As was shown analytically in Ref-j^l, nonlinear lattices of the FPU type have an impor- 
tant peculiarity. Specifically, the unperturbed motion is linear, therefore, the KAM theory 
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mally cannot be applied. Indeed, one of the conditions for applicability of 
the KAM theory is that the unperturbed frequency of oscillations has to be dependent on 
energy. For this reason, there is no rigorous ground to expect that under a very weak per- 
turbation the motion of the FPU model remains stable. Indeed, the analysis of Ref. [89( 
has shown that the resonance overlap can occur for arbitrarily small but nonzero values of 
the parameters a and (3. Therefore, stochastic motion does not disappear in the limit of 
vanishing perturbation, although the degree of stochasticity may be very small. This situa- 
tion, known as "nearly linear oscillations" arises in many practical applications and requires 
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specific methods of study (9^]. 

Another effect which is important in view of discussion about the behavior of the FPU 
models for a weak perturbation, is the influence of nonlinear resonances that appear in high 
orders of perturbation theory. Indeed, the overlap criterion [f| is based on consideration of 
the first-order resonances only. On the other hand, second-order resonances may create a 
wide resonance net and lead to strong equipartition which can occur at much larger time 
scales. This problem was briefly discussed in Refs.jgl Q], and further numerical studies 
have confirmed the point that apart from strong chaos (with a fast equipartition and large 
values of the Lyapunov exponents), one can speak about weak chaos. This weak chaos is 
characterized by a much weaker equipartition of energy, smaller values of LE, and a different 
kind of scaling of the LE. Recently, the role of these high-order resonances was examined 
in Ref.jsjj], together with a detailed analysis of peculiarities related to the nearly-linear 
character of the oscillations. The a— model was examined and the stochasticity threshold 
has been found for both strong and weak interactions. 

The concept of weak chaos poses an important question of whether the threshold of 
stability for an infinite time scale vanishes in the thermodynamical limit, N — > 00. To 
shed some light in this problem, in Ref.j^] numerical simulations have been made for the 
(3— model (with periodic boundary conditions) paying main attention to the long-time be- 
havior. Measuring the spectral entropy ([22)1 . it was found that the result strongly depends 
on whether initially one mode, Ak = 1 or a group of modes with Ak 3> 1, is excited (with 
both k and Ak proportional to N). In the first case, the data indicate the existence of a 
threshold (3E/N ~ const. In the second case the threshold vanishes as (3E/N ~ 1/iV 7 with 
7 ph 1. As one can see, these results are in contradiction with the estimate (fTSj) obtained 
for resonance overlap. The existence of finite times of relaxation to the equipartition has 
also been confirmed in Ref . [93] where the strong influence of initial transient times has been 
detected, after which generic behavior emerges with no dependence on initial conditions. 

One specific question widely discussed in the literature is how fast the relaxation is to a 
steady-state distribution of energy among normal modes. An extensive study of this problem 



was performed in Ref. 



see, also, Ref.|94|). The authors made an attempt to relate the 



data for time dependence of the spectral entropy t](t), to the estimates of Refs. |95l l96l . IqT | 
for the rate of the Arnold diffusion. In 1964 Arnold has proved 9^ that many- dimensional 
nonlinear systems are, in general, globally unstable due to a very peculiar diffusion (Arnold 
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diffusion), for details see, e.g. Refs.Il2l Q). Loosely speaking, this diffusion occurs (below 
the border of resonance overlap) for initial conditions inside the narrow stochastic layers 
which surround any nonlinear resonance. Due to an everywhere dense set of resonances (of 
different orders), starting from a point inside this Arnold web, the trajectory diffuses over 
the web along these resonances. Although Arnold diffusion is extremely weak (exponentially 
small in the perturbation parameter), the motion is unbounded in the phase space of the 
system. There is a widespread belief (although still there are no reliable numerical results) 
that Arnold diffusion is responsible for a weak instability in FPU lattices. By fitting data to 
Nekhoroshev's expressions, in Ref.jsi| the following empirical dependence has been obtained: 

rj(t) = exp [- (t/rH (27) 

for t < t r , and 

ri(t) = r/oo = exp [- {tr/tY} (28) 

for t > tr, with the numerical estimate of 0.3 < v < 0.5. As for the relaxation time 
tr, it is argued that it is proportional to the energy density, r R ~ e, and therefore, to the 
number of particles N for fixed total energy. The nonzero value of is discussed in Ref. 
Among others, the most realistic explanation of this result is that it is due to fluctuations 
of the energies Ek, which are not taken into account in the normalization factor for 77. 
One of the important conclusions drawn in Refs. is that the equipartition of energy 

is always reached. This supports the expectation of non-existence of a minimal critical 
value of nonlinearity for the stochasticity. Another conclusion is that the critical value of 
perturbation e c which marks the transition from weak to strong stochasticity, corresponds 
to the overlap condition (fTBj) . 



The region of weak stochasticity was closely examined in Ref. |10f| . For initial conditions 
the low modes were excited in the (3— model with zero boundary conditions. As was shown 
numerically, for low initial energy the distribution of mode energies after large time follows 
an exponential decrease with increasing the mode number. This exponential dependence was 
explained theoretically, by finding an approximate solution of equations of motion written 
in a form similar to that of a nonlinear Schrodinger equation. It was also shown numerically 
that a single nonlinearity parameter, 

R = Qtt- 2 (3E 1 (N + 1), (29) 
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governs the local interactions between the low k N modes. For R ^> 1 there is a critical 
energy E c « 2.8 (for /3 = 0.1) for which strong diffusion in k— space leads to equipartition. 
Below this border the diffusion leads to an exponential distribution of mode energies. It 
is noted that R being sufficiently large corre spon ds to the overlap criterion of the onset of 



widespread stochasticity. In further studies 



10l| | the energy transitions and different time 



assifi ed and discussed, as a function of energy. These results have been 



102] to FPU- type lattices with two types of masses randomly distributed 



scales have been c 
generalized in Ref. 
along the chain. 

The dependence on initial conditions is another important question. One can classify the 
following cases: (i) single mode excitation with small k; (ii) a group of low frequency modes 
with Ak and k proportional to iV (thermodynamic limit); (iii) single mode excitation with 
k « N/2 (narrow packet approximation), and (iiii) high frequency excitation with k « N. 
As one can see, other important cases are missed in this picture, thus showing how difficult, 
in general, the problem is of stati stica l relaxation in nonlinear lattices. A comparison of 
cases (i) and (ii) was done in Ref. |l0l| . It was found that transient times to equipartition 
in case (ii) are proportional to y/~N, in comparison with single mode excitations (i) where 
this time does not exist or is not important. The transient times are characterized by non- 
universal dynamical characteristics, in contrast with larger times for which one can find a 
good scaling dep endence on the energy density E/N. Concerning the case (iiii), one can 
refer to Ref. [103J where the energy equipartition for initially excited high frequency modes 
was studied. As for the case (iii), the important findings are analytical ones, with much 
attention given to the instability conditions. 



VIII. NSE AND BEC 



As was found in Ref. |69j, one of remarkable properties of the (3— model (J2J) is its direct 
relevance to the nonlinear Schrodinger equation (NSE). Indeed, let us write the Hamiltonian 
corresponding to the equations of motion (J2J), 



N 



H = J2 



n=l 



V 2 1 

2 2 v 



(30) 
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In the following, we consider periodic boundary conditions, xq = xn- In analogy with the 



quantum mechanics, we use the canonical variables a k and a 



where 



ak = ^ = k {Vk ~ tUJkQ * k) (31) 

i JV i N 

^ = -Jjr T £ Vne~ % — ■ & = -7= E ^ e ~ ( 32 ) 

V^V n =l V n =l 

with = 2 sin ^ as the frequency of the k— th linear mode, = 1, 2, iV. Assu ming that 
the initial packet in k— space is narrow and centered at k « AT/2, it can be shown 69] that 
the Hamiltonian takes the form, 

H = ^ UJka*k a k + ^ ^kik2k 3 k i ^k 1 a k2 a k-3 a k^k 1 +k2-k- A -k i + 0(1). (33) 

fc=l kik2ksk4 



Here the term = V + W (gi + <? 2 + <?3 + 94) with V =_|Ll sm7r Ar) and W = 

l^sin^?^ <C Vo describes the (resonant) four- wave interaction [104| (with q = k — k <C 
k ). All other (non-resonant) terms of 0(1) can be neglected in this approximation. By 
expanding u>k at the point ko, one can obtain, Uk ~ ujk + Aq — Qq 2 . The parameters A and 
Q are A = ^ cos ^ « (^/A^) 3 and = ) 2 sin ^ « (tt/A^) 2 . As a result, the equations 
of motion ihk = dH / da* k can be written as 

iA q = -Qq 2 A q + V J2 K A ^A+n-n-j s (34) 

9l,?2,93 

where A q = exp (i(^fc + Ag)t) a g +fc . As one can see, these equations describe a nonlinear 
chain of interacting oscillators. Using the transformation <&(0,t) = J2 q A g (t) exp(iq6) = 
&(9 + 2tt, t), we obtain the NLS equation, 

i _ = n 7jF + V.|»r». (35) 

This classical equation is well known in the physics of interacting particles and widely 
discussed in many applications. It is a particular case of the Gross-Pitaevskii (GP) equation 
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1061 ] which attracts much attention in connection with Bose-Einstein condensation 
(BEC). Using the mean-field approximation, this equation describes the evolution of the 
condensate wave function, and, in essence, is of a semi-classical nature. Th e important 
peculiarity of the GP-equation is its complete integrability (see, e.g., Ref. jl07 |). This is of 
special interest from the point of view of the statistical properties of the condensate. As was 
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confirmed numerically 69], in the FPU (3— model good conservation of first three integrals 
of motion, analytically derived for the GP-equation, can be observed provided the initial 
packet is centered at k = N/2. 

The model of type (J33|) was recently examined using close numerical investigations 

jl08^ . Specifically, the dynamics of the condensate in one-dimensional geometry was as- 
sumed to be governed by the following Hamiltonian in the action-angle representation, 
A n = A/Inexp (i6 n ): 

# = 5> n / n +f £ K 1 n 2 n3n 4 (/n 1 /n 2 /n 3 /n 4 ) 1/2 e-<^ +e --^-^). (36) 

Here htu n = n 2 7r 2 h 2 /2mL 2 labels the single-particle energy levels, m is the mass of particles, 
L is the length of the one- dimensional box, and V nin2n3ri4 corresponds to the matrix elements 
Jo <t>m 4>no4>n -,4>nAdz of the interaction, with (p n as the normalized modes of the box (see details 
in Ref.[l08|). 

In the study of the model ()36|). methods developed for classical nonlinear lattices have 
been extensively used. Of particular interest was the time-dependence of the normalized 
spectral entropy since it can be used to distinguish between regular and irregul ar d^ 
namics of the condensate. Note that due to the zero boundary conditions used in Ref. 
the dynamics may be non-integrable in contrast to the GP-equation. For a weak interac- 
tion the dynamics was numerically found to be quasi-periodic for generic initial conditions, 
which may be compared with the recurrent motion in the FPU models. On the other hand, 
with increased interaction strength, the dynamics appears to reveal chaotic properties. This 
effect was accompanied by a strong decrease of the spectral entropy f](t), and by irreversible 
dynamics. The main result is that starting from conditions in which the low modes were 
initially excited, the energy diffusively spreads over modes with higher frequencies. The 
stochasticity was found to be stronger for low excited modes, in contrast to the FPU model. 
This fact may be explained by the different kinds of the unperturbed frequency spectrum in 
these two models. 

The above results show that many of properties of classical nonlinear lattices, especially, 
non-integrable models of the FPU type, can be discussed in a more general context, namely, 
in the context of the physics of interacting quantum particles. Although the latter subject 
is not new, the problem of the onset of many-body chaos in quantum systems of interacting 
Fermi and Bose - particles has recently attracted much attention due to its important phys- 
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rdy- 



Ref. 



Jifl 



ical appl i catio ns (see, e.g. Ref. [lQ9|). In this respect, an interesting study was reported in 



111] where the quantum analog of Eq.(p£i)) was introduced, 



J2J3J4 



J+J2-J3-J4 



(37) 



Aj, A\ 



Here 
model. 

Analytical studies 



5jk and q = h(3 cot(7r/2iV)/32iV, with u and Vq defined as in the classical 
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1 111 ] of the dynamical instability of motion gave quite unexpected 



results. As was already discussed in Section VI, in the narrow packet approximation the 
classical dynamics of the FPU model displays an instability above the critical value given 
by Eq.(|25)). In contrast to this result supported by other studies, in the quantum model the 
instability occurs for any weak perturbation. On the other hand, the rate of this instability 



je explai ned 
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^y 



turns out to be slower than for the classical model. If the latter effect could 
quantum suppression of classical instability (as in the Kicked Rotor model 
the absence of the instability threshold is a somewhat new effect. One can suggest that this 
effect is due to quantum tunneling in the effective potential, and further studies appear to 
be very important. 



The above approach has been recently developed and applied in Ref. 115] to the problem 
of collapse in the Bose-Einstein condensation with an attractive potential. As is known, 
the collapse dynamics cannot be described b y the Gross-Pitaevskii equation, therefore, new 



ideas are important. As was shown in Ref. 



115J, near the instability threshold quantum 



effects turn out to play an important role and have to be taken into account. Specifically, 
when comparing the GP and quantum dynamical growth of the unstable modes, the absence 
of the instability threshold was confirmed for the quantum model, in contrast with the GP 
equation. This means that the quantum solution is always unstable, and eventually collapses 
in a finite time. The difference between the GP and the quantum model is important when 
approaching the classical (GP) critical limit of the instability. This effect may be compared 
with an exponentially fast spread of wave packets in quantum systems which are chaotic in 
the classical limit. As was shown in Refs. jlld Ill7^ . for classically chaotic systems quantum 
effects are extremely str ong and reveal themselves on a very short (logarithmic) time scale. 
Recent rigorous results |ll8j | on the quantum instability of averages near stationary points 
seem to have a direct relevance to the more general problem of quantum corrections in the 
region of chaos. 
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It would be interesting to study the influence of terms, neglected in the narrow packet 
approximation, when obtaining the NLS equation (|35|). It is easy to show that by keeping 
the next terms in the expansion of u k and Vk u k 2 ,k 3 ,k4 about the point k , one can derive the 
following equation: 

i _ = n _ + V!l |,|^ + ix __ Wo |,|._, (38) 

which describes the dynamics of the FPU model more correctly than the Gross-Pitaevskii 
equation (with \ — ^d 3 Uk/dk 3 at k — ko). The additional terms in the above equation may 
be important for describing the evolution of wave packets for large time scales, and for the 
breakdown of integrability. 

One recent attempt to consider the dynamics of the Bose-Einstein condensate by m aking 
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use of an approach developed in the field of quantum chaos, was performed in Ref. 
Using numerical simulations, the authors focus on the dynamics of the Bose-Einstein con- 
densate on a ring, which is described by the quantum Hamiltonian, 

H = J2e k h k + ^- a{ala p a r 5 k+q - p - r . (39) 

k k,q,p,r 

Here h k = a\ak is the occupation number operator, and a s are the creation-annihilation 
operators, and e k = 47r 2 k 2 /L 2 . As one can see, this Hamiltonian can be considered to be 
the quantum version of the classical Hamiltonian (|33J) describing the evolution of the FPU 
model in the narrow packet approximation. The behavior of this system is governed by 
only one parameter n/g, where n is the particle density on a ring of length L, and g is the 
strength of the interaction between bosons, determined by the interatomic scattering length. 

As is known, for weakly interacting particles, n/g — > oo, the mean-field approximation 
gives the correct description of the dynamics. In the other limit of strongly interacting 
particles, known as the To nks- Girardeau regime, n/g —>■ 0, the density of interacting bosons 
becomes identical to that of non-interacting fermions. The transition between these two 
regimes is known to corresp ond, approximately, to n/g ~ 1. 

The main interest in Ref. [HS| was to observe and quantify the degree of irregularity in the 
dynamics of the condensate. Specifically, the situation in which all bosons initially occupy 
the single-particle level with the angular momentum k = has been explored, with a fu rther 



analysis of the evolution of the condensate in time. The main result of the study in Ref. 119 1 



was that with an increase of the interaction strength g, regular (quasi-periodic) dynamics 
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alternates with irregular behavior of the observable quantities. This transition was found 
numerically to occur at the transition from the mean-field to the Tonks-Girardeau regimes. 
Given the clear evidence of the efficiency of the proposed approach, these results open the 
door for further studies of the condensate dynamics from the viewpoint of many-body chaos. 



IX. CONCLUDING REMARKS 



Due to space limitations, many recent studies are not discussed here, although they are 
relevant to the FPU model. In the first line, one should mention an increasing interest in 
the existence of localized n onlin ear oscillations (breathers) emerging in nonlinear lattices 
(for a review, see, e.g. Ref. 120|). For some time, their existence in the FPU model was 
questionable, mainly due to the fact that the main int erest initially was related to low- 
frequency excitations. As was shown in the early paper 1211 ] . localized optical excitations 
(high-frequency modes) can be observed in the FPU model with alternative masses. This 
was the first indication that by exciting the highest modes in the FPU lattice, a new kind 
of solution with special str uctu re emerges (the most recent results on the diatomic FPU 



model are reported in Ref. 



122j). Although self-localized solit ons in an harmonic lattices 



without impurities were predicted quite a long time ago in Refs. 



123 



124j , only re cent 



existence of breathers in FPU lattices has been proved rigorously (see, e. g., Refs. 



y the 



125j). In 



where it 



connection with the FPU problem one should mention the results of Refs. 
was shown that breathers can be responsible for the slow relaxation of initially thermalized 
nonlinear lattices. This effect seems to be relevant to the long-term regular dynamics in the 
FPU models, which is alternated by a strong energy sharing betwee n lin ear modes. To date, 
many studies of breathers (including chaotic breathers, see in Ref. jl28^ ) in the FPU model 
have been carried out, and we hope that the reader can find in this issue more information 
on the subject. 

Coming back to the original question about the ergodicity and thermalization in the 
FPU model, one should conclude that some of problems still remain open. In particular, the 
existence of the threshold for weak chaos in the thermodynamical limit N —>■ 00 is still under 
study by many researchers. As is clear from the above discussions, the main difficulty, apart 
from the numerical one, is the strong dependence of the results on the model parameters. 
The behavior of the model depends strongly on whether low- or high-frequency modes are 
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initially excited. Also, the number of excited modes seems to be important for the dynamics, 
as is indicated in previous studies. Last, but not least, is the fact of a difference between 
the alpha— and beta— models. Therefore, future studies are desirable, both analytical and 
numerical. 



One should mention the direct relevance o 
Einstein condensation. As was shown in Ref. 



t 



re FPU model to the models of the Bose- 



69(, the narrow packet approximation in the 



(3— model leads to the Gross-Pitaevskii equation with an attractive potential. Thus, the 
instability of highest modes in periodic FPU lattices corresponds to that of the Bose-Einstein 
condensate. This fact is important for further studies of instabilities both in the nonlinear 
classical lattices and in quantum models of the Bose-Einstein condensate. Note that now 
it is possible to control experimentally the sign of the intera ction between bosons, and to 
observe the collapse of the condensate (see, for example, Ref.|l29|). 

Another new direction is the study of dynamics of quantum models of interacting Bose 
particles. In particular, direct quantization of classical nonlinear chains related to the Gross- 
Pitaevskii equation sho ws a close analogy for the dynamics in classical and quantum models. 
Recent numerical data jl 1 o| for the transition between the mean field and Tonks- Girardeau 
regimes have revealed the onset of irregular motion of the condensate. This type of transition 
is known to occur in quantum models of interacting particles which are chaotic in the classical 
limit. Therefore, one can expect that methods well developed in the theory of quantum 
chaos, may give new insight on the dynamics of interacting bosons in the condensates. 

The above practical problems are part of the more general problem of quantum-classical 
correspondence for systems with irregular behavior. From this point of view, the FPU model 
is o f par ticular interest and can be considered as an important example. As discussed in 
Ref. 130], there are two mechanisms which are responsible for the appearance of statistical 
behavior of dynamical (deterministic) systems. The first mechanism is the thermodynamic 
limit with iV — > oo, which is well known since the early days of statistical mechanics. The 
important point here is that this limit has nothing to do with chaos, it is based on the 
ergodicity of motion only, which is known to be the weakest statistical property. As we 
already discussed, perfect statistical and thermodynamical properties are known to emerge 
even in completely integrable systems such as the Toda-lattice. Although there are initial 
conditions which correspond to solitons, the measure of these specific conditions is extremely 
small, and surely can be neglected practically. Therefore, the role of additional terms that 
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break the integrability becomes important when the number N is finite. The expectation 
of FPU was based on their belief that N = 32, 64 is large enough in order to observe the 
equipartition in the presence of small nonlinearity. 

The second mechanism for the onset of statistical properties in dynamical systems is, in 
principle, different. It is based on a local instability of motion for generic initial conditions 
in the phase space of the system. With this mechanism the ergodicity is not important 
provided the total measure of initial conditions with regular motion is very small, although 
it can be finite. Due to this local instability (with reflecting boundaries in phase space), 
this motion reveals clear mixing properties, leading to strong sensitivity of the motion to 
initial conditions. As a result, an apparent irreversibility of motion emerges since any weak 
external perturbation gives rise to non-recurrence of the initial conditions. The important 
point for this scenario is that the dynamical chaos can emerge in systems with few degrees of 
freedom, in contrast to the first (thermodynamic) mechanism. It is important to stress that, 
although these two mechanisms are different, the common feature is that in both cases the 
time dependence of the observables can be desc ribed by an infinite number of statistically 
independent frequencies (see details in R.ef . jl3n| ) . 

Turning to quantum systems, the origin of "quantum chaos" in dynamical systems is 
based on the first mechanism, with no relevance to the local instability of trajectories. Specif- 
ically, irregular behavior of a system emerges when an initial wave packet consi sts o f many 
exact chaotic eigenstates with statistically independent frequencies (see, e.g. Ref. |l3l| ). This 
concept is very important for establishing the conditions for the onset of chaos in quantum 
systems and in quantifying their irregular properties. Therefore, the understanding of phys- 
ical effects found in the FPU model, as well as the use of tools developed for identifying 
these effects, may be useful for the study of quantum systems with complex behavior. 

As a consequence of their research on the FPU model, in this brief review the authors 
have tried to summarize the main ideas, tools and results related to the FPU paradox, 
after 50 years of the celebrated paper. For one of the authors (FMI), the FPU problem 
was his initial scientific PhD research. For GPB, the relation discovered between the FPU 
model and nonlinear Schrodinger equation has contributed to his interests in Bose-Einstein 
condensation problems. And for both of the authors, the preparation of this manuscript 
provided a welcome opportunity to learn much more about new trends in the physics of 
nonlinear phenomena. 
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